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ABSTRACT 

This  paper  examines  the  limitations  of  using  B -spline  representation  as  an  analysis  tool  by 
comparing  its  geometry  with  the  nonlinear  finite  element  absolute  nodal  coordinate  formulation 
(ANCF)  geometry.  It  is  shown  that  while  both  B -spline  and  ANCF  geometries  can  be  used  to 
model  non- structural  discontinuities  using  linear  connectivity  conditions,  there  are  fundamental 
differences  between  B-spline  and  ANCF  geometries.  First,  while  B-spline  geometry  can  always 
be  converted  to  ANCF  geometry,  the  converse  is  not  true;  that  is,  ANCF  geometry  cannot  always 
be  converted  to  B-spline  geometry.  Second,  because  of  the  rigid  structure  of  the  B-spline 
recurrence  formula,  there  are  restrictions  on  the  order  of  the  parameters  and  basis  functions  used 
in  the  polynomial  interpolation;  this  in  turn  can  lead  to  models  that  have  significantly  larger 
number  of  degrees  of  freedom  as  compared  to  those  obtained  using  ANCF  geometry.  Third,  in 
addition  to  the  known  fact  that  B-spline  does  not  allow  for  straight  forward  modeling  of  T- 
junctions,  B-spline  representation  cannot  be  used  in  a  straight  forward  manner  to  model 
structured  discontinuities.  It  is  shown  in  this  investigation  that  ANCF  geometric  description  can 
be  used  to  develop  new  spatial  chain  models  governed  by  linear  connectivity  conditions  which 
can  be  applied  at  a  preprocessing  stage  allowing  for  an  efficient  elimination  of  the  dependent 
variables.  The  modes  of  the  deformations  at  the  definition  points  of  the  joints  that  allow  for  rigid 
body  rotations  between  ANCF  finite  elements  are  discussed.  The  use  of  the  linear  connectivity 
conditions  with  ANCF  spatial  finite  elements  leads  to  a  constant  inertia  matrix  and  zero  Coriolis 
and  centrifugal  forces.  The  fully  parameterized  structural  ANCF  finite  elements  used  in  this 
study  allow  for  the  deformation  of  the  cross  section  and  capture  the  coupling  between  this 
deformation  and  the  stretch  and  bending.  A  new  chain  model  that  employs  different  degrees  of 
continuity  for  different  coordinates  at  the  joint  definition  points  is  developed  in  this 
investigation.  In  the  case  of  cubic  polynomial  approximation,  C1  continuity  conditions  are  used 
for  the  coordinate  line  along  the  joint  axis;  while  C°  continuity  conditions  are  used  for  the  other 
coordinate  lines.  This  allows  for  having  arbitrary  large  rigid  body  rotation  about  the  axis  of  the 
joint  that  connects  two  flexible  links.  Numerical  examples  are  presented  in  order  to  demonstrate 
the  use  of  the  formulations  developed  in  this  paper. 

Keywords:  Geometric  discontinuities;  Finite  element;  Multibody  systems;  B-spline;  NURBS. 
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1.  INTRODUCTION 

The  geometry  description  used  in  many  of  the  existing  finite  element  (FE)  formulations  cannot 
be  exactly  converted  to  the  geometry  developed  by  computational  geometry  methods  such  as  B- 
spline  and  NURBS  (Non-Uniform  Rational  B-Splines)  representations.  This  fact  has  motivated 
researchers  in  the  mechanics  community  to  adopt  the  methods  of  computational  geometry  as 
analysis  tools  instead  of  using  conventional  FE  formulations.  While  the  methods  of 
computational  geometry,  such  as  B-spline,  have  several  desirable  analysis  features;  these 
methods  have  serious  limitations  when  used  as  analysis  tools.  The  B-spline  recurrence  formula 
and  the  rigid  definition  of  the  knot  vector  make  B-spline  less  attractive  as  compared  to  the 
absolute  nodal  coordinate  formulation  (ANCF)  geometry  description.  While  B-spline  geometry 
can  always  be  converted  exactly  to  ANCF  geometry  (Piegl  and  Tiller,  1997,  Sanborn  and 
Shabana,  2009;  Fan  and  Shabana,  2010),  the  converse  is  not  always  true.  ANCF  geometry  does 
not  restrict  the  order  of  the  parameters  or  the  number  of  basis  functions  used  in  the  interpolating 
polynomials  (Dmitrochenko.and  Pogorelov,  2003;  Dufva  et  al.,  2005;  Garcia- Vallejo  et  al., 
2008;  Kerkkanen  et  al.,  2006;  Schwab  and  Meijaard,  2010;  Tian  et  al.,  2009,  2010;  Yoo  et  al., 
2004;  Yakoub  and  Shabana,  2001;  Shabana  and  Mikkola,  2003;  Abbas  et  al.,  2010).  This 
advantage,  as  will  be  demonstrated  in  this  paper,  allows  for  developing  finite  elements  with  less 
number  of  degrees  of  freedom  as  compared  to  those  developed  using  the  B-spline  geometry. 
Another  fundamental  difference  between  B-spline  and  ANCF  geometric  descriptions  lies  in 
modeling  discontinuities.  As  previously  explained  by  the  authors  using  simple  planar  examples, 
there  are  two  types  of  discontinuities  when  chain  systems  are  considered  (Flamed  et  al.,  2011). 
The  first  is  structural  discontinuity  which  does  not  allow  rigid  body  displacement  between  two 
elements  connected  at  the  joint  definition  point.  This  joint  allows  only  for  deformation  degrees 


of  freedom.  The  second  type  of  discontinuity  is  called  non-structurcd  discontinuity  which  allows 
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for  rigid  body  displacement  at  the  joint  definition  point.  Figure  1  shows  a  chain  which  has  a 
structural  discontinuity  at  point  C  and  non-structural  discontinuity  at  point  O  .  At  the  junction  at 
C,  only  deformation  degrees  of  freedom  are  allowed,  while  at  point  O,  relative  rigid  body 
rotation  is  permitted.  Nonetheless,  the  degree  of  continuity  at  both  points  is  C° .  B-spline  can  be 
used  as  an  analysis  tool  to  describe  the  non-structural  C°  discontinuity  at  point  O ,  but  because 
of  its  rigid  recurrence  structure  and  the  definition  of  its  knot  vector  and  knot  multiplicity,  B- 
spline  cannot  be  used  in  the  motion  analysis  of  structural  C°  continuity  at  point  C  since  B- 
spline  C°  description  leads  to  a  rigid  body  mode;  that  is,  the  elimination  of  one  control  point  by 
reducing  the  knot  multiplicity  by  one  is  not  sufficient  for  eliminating  the  modes  of  rigid  body 
rotations  between  two  B-spline  segments.  ANCF  geometry,  on  the  other  hand,  can  be  used  in  the 
analysis  of  both  structural  and  non-structural  discontinuities  (Hamed  et  al.,  201 1). 

One  of  the  important  multibody  system  (MBS)  applications  that  can  be  used  to  shed  light 
on  the  fundamental  differences  between  B-spline  and  ANCF  geometries  are  chain  applications. 
Chains  are  highly  nonlinear  systems  that  are  subjected  to  repeated  impulsive  forces  during  their 
functional  use.  Geometric  nonlinearities  are  the  result  of  the  large  relative  displacements 
between  the  chain  links.  The  repeated  impulsive  forces  as  the  result  of  the  chain  link  contact  with 
the  rollers  and  other  system  components  introduce  high  frequencies  to  the  nonlinear  chain 
dynamic  model.  For  these  reasons,  the  nonlinear  dynamic  analysis  of  chain  systems  represents 
one  of  the  most  challenging  computational  problems.  In  fact  the  simplest  rigid-link  chains  are 
highly  nonlinear  because  of  the  large  relative  rotations  between  the  chain  links  (Roberson  and 
Schwertassek,  1988;  Schiehlen,  1997).  In  the  case  of  rigid-link  chains,  the  geometric 
nonlinearities  that  result  from  these  rotations  lead  to  highly  nonlinear  chain  inertia  forces  that 
include  the  quadratic  velocity  Coriolis  and  centrifugal  forces.  Furthermore,  the  contact  between 
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the  rigid-chain  links  and  the  sprockets  and  rollers  that  may  exist  in  the  system  is  often  described 
using  compliant  force  models;  leading  to  high  frequencies  that  require  the  use  of  very  small  time 
integration  step  in  order  to  accurately  capture  the  changes  in  the  velocities,  accelerations,  and 
forces.  For  these  reasons,  efficient  and  accurate  modeling  of  chain  systems  remains  a  challenging 
MBS  computational  problem  even  in  the  simpler  case  in  which  the  chain  links  are  assumed  to  be 
rigid. 

Because  of  the  geometric  nonlinearities  and  the  high  frequencies,  it  is  important  to  use  an 
efficient  solution  algorithm  if  the  flexibility  of  the  chain  links  is  considered.  Flexible-link  chains 
require  the  use  of  significantly  larger  number  of  degrees  of  freedom  in  order  to  capture  the  link 
deformation  modes.  Some  of  these  deformation  modes  may  also  introduce  high  frequencies  in 
addition  to  the  high  frequencies  resulting  from  the  contact  between  the  chain  links  and  rollers  as 
well  as  other  components  in  the  system.  It  is  also  important  in  some  applications  to  capture 
certain  coupled  deformation  modes  that  cannot  be  captured  using  conventional  structural  finite 
elements  such  as  beams  and  plates  that  are  based  on  simplified  kinematic  assumptions.  For 
example,  in  tracked  vehicle  applications,  the  chain  links  are  subjected  to  significant  tensile  and 
compressive  forces.  The  coupling  between  the  deformation  of  the  link  cross  section  and  other 
modes  of  deformation  can  be  significant  and  must  be  taken  into  account  in  order  to  develop  a 
more  realistic  model.  It  is,  therefore,  important  to  employ  FE  formulation  that  captures  the  effect 
of  these  coupled  deformation  modes  and  allows  for  an  efficient  MBS  implementation.  It  was 
shown  in  the  literature  using  planar  examples  that  ANCF  finite  elements  can  be  used  to  develop 
new  FE  meshes  for  chain  applications  (Garcia  de  Valljo  et  al,  2003;  Hamed  et  al,  2011).  In  these 
FE  meshes,  the  flexible-link  pin  joints  are  defined  using  linear  connectivity  conditions  despite 
the  large  relative  rotation  allowed  between  the  chain  links.  This  leads  to  an  efficient  elimination 
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of  the  dependent  variables  at  a  preprocessing  stage.  Furthermore,  the  use  of  the  linear 
connectivity  conditions  with  ANCF  finite  elements  leads  to  a  constant  chain  inertia  matrix  and 
zero  Coriolis  and  centrifugal  forces.  It  is  one  of  the  main  objectives  of  this  investigation  to 
demonstrate  for  the  first  time  that  a  three-dimensional  flexible-link  chain  model  that  is  based  on 
linear  connectivity  conditions  and  has  a  constant  mass  matrix  and  zero  Coriolis  and  centrifugal 
forces  can  be  developed  using  spatial  fully  parameterized  ANCF  finite  elements. 

This  paper  focuses  on  two  fundamental  issues  that  summarize  its  main  contributions. 
First,  the  fundamental  differences  between  B -spline  and  ANCF  geometries  are  demonstrated.  It 
is  shown  that  while  B-spline  geometry  can  always  be  converted  to  ANCF  geometry,  the  converse 
is  not  true  because  of  the  rigid  B-spline  recurrence  structure.  It  is  also  shown  that  B-spline 
representation  can  be  used  only  in  the  analysis  of  one  type  of  C°  continuity  referred  to  in  this 
paper  as  non-structural  discontinuity;  while  such  a  B-spline  representation  cannot  be  used  in  the 
analysis  of  another  type  of  C°  continuity  referred  to  as  structural  discontinuity.  It  is  shown  that 
ANCF  finite  elements  which  have  degrees  of  freedom  less  than  their  B-spline  counterparts  can 
be  developed  since  ANCF  does  not  have  specific  requirements  on  the  order  of  the 
parameterization  variables  or  the  number  of  basis  functions  used  in  the  interpolating 
polynomials. 

The  second  main  contribution  of  this  paper  is  to  develop  a  new  three-dimensional 
flexible-link  chain  model  using  fully  parameterized  ANCF  finite  elements.  This  chain  model  is 
based  on  a  new  FE  mesh  defined  using  linear  connectivity  conditions.  The  FE  mesh  allows  for 
relative  rigid  body  rotations  between  its  elements  and  has  a  constant  inertia  matrix  and  zero 
Coriolis  and  centrifugal  forces.  In  order  to  develop  the  new  flexible-link  chain  model  presented 
in  this  paper,  a  new  pin  joint  model  is  introduced.  At  the  joint  definition  point,  different  degrees 
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of  continuity  are  used  with  different  parameters;  leading  to  some  strain  components  to  be 
continuous  while  the  others  are  discontinuous.  The  modes  of  deformation  at  the  joint  definition 
points  are  discussed  in  order  to  shed  light  on  the  nature  of  the  new  joint  and  kinematic 
constraints  developed  in  this  paper.  Numerical  results  are  presented  in  order  to  demonstrate  the 
use  of  ANCF  finite  elements  in  developing  the  new  flexible-link  chain  model.  The  limitations 
identified  in  this  paper  when  B-spline  geometry  is  used  as  analysis  tool  suggest  the  use  of  the 
integration  of  computer  aided  design  and  analysis  (I-CAD-A)  approach  (Sanborn  and  Shabana, 
2009;  Lan  and  Shabana,  2010).  In  the  I-CAD-A  approach,  a  constant  transformation  can  be 
developed  to  convert  CAD  geometry  to  FE  mesh. 


2.  B-SPLINE  SURFACES 

B-spline  surfaces  are  defined  using  the  product  of  B-spline  base  functions,  two  parameters,  and 
two  knot  vectors.  B-spline  surfaces  can  be  defined  in  the  following  parametric  form  (Piegl  and 
Tiller,  1997): 

n  m 

«'= 0  7=0 

where  u  and  v  are  the  parameters;  Nj  p{u)  and  Njq(v )  are  B-spline  basis  functions  of  degree 
p  and  q ,  respectively;  and  P  .  are  a  set  of  bidirectional  net  of  control  points.  The  B-spline 
basis  functions  Nj  p(u )  are  defined  as 
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N,  „(«)  =  ■ 


1  if  ui  <u<  u 
0  otherwise 


i+i 


Nu(u)  =  — — —Ni  hl(u)  +  Ui+i+l  -  N. 


Ui+j  —  U; 


M,+;+I  ~UM 


i+hj- 


1  («) 


(2) 


where  i  =  0, 1,2,...,/;  +  p  + 1  are  called  the  knots',  and  ut  <  uM .  The  vector 


U  =  {u0  ux  un+  x }  is  called  the  knot  vector.  Similar  definitions  can  be  introduced  for 


Nj  (y)  with  another  knot  vector  V  =  {v0  v,  •••  vm+q+x}.  Note  that  the  orders  of  the 

polynomials  in  the  u  and  v  directions  can  be  different;  for  example,  a  cubic  interpolation  can  be 
used  along  u  while  a  linear  interpolation  can  be  used  along  v .  As  in  the  case  of  B -spline  curves, 
the  knots  of  B-spline  surfaces  do  not  have  to  be  distinct;  distinct  knots  are  called  breakpoints  and 
define  surface  segments  with  non-zero  dimensions.  The  number  of  the  non-distinct  knots  in  U 
and  V  at  a  point  is  referred  to  as  the  knot  multiplicity  associated,  respectively,  with  the 
parameters  u  and  v  at  this  point.  At  a  given  breakpoint,  the  multiplicity  associated  with  u  can 
be  different  from  the  multiplicity  associated  with  v  ;  allowing  for  different  degrees  of  continuity 
for  the  derivatives  with  respect  to  u  and  v.  For  cubic  Nip  (p  =  3),  C°,C',  or  C2  conditions 
correspond,  respectively,  to  knot  multiplicity  of  three,  two,  and  one;  while  in  the  case  of  linear 
interpolation  of  AF  ,  the  highest  continuity  degree  that  can  be  demanded  is  continuity  of  the 

gradients.  When  zero  multiplicity  is  used  at  a  breakpoint,  the  segments  blend  together  at  this 
point. 

In  B-spline  surface  representation,  there  is  a  relationship  between  the  polynomial  degree, 
the  number  of  knots,  and  the  number  of  control  points.  This  relationship  must  be  fully 
understood  if  B-spline  geometry  will  be  used  as  an  analysis  tool.  If  r  + 1  is  the  number  of  knots 
in  U  and  5  +  1  is  the  number  of  knots  in  V  ,  then  in  B-spline  geometry,  one  must  have 
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r  =  n  +  p  + 1,  s  =  m  +  q  + 1  (3) 

These  formulas  imply  that,  for  a  given  polynomial  order,  if  the  number  of  knots  decreases,  the 
number  of  control  points  (degrees  of  freedom  used  in  the  analysis)  must  also  decrease.  A 
decrease  in  the  knot  multiplicity  by  one  is  equivalent  to  eliminating  one  control  point.  This  can 
also  be  equivalent  to  increasing  the  degree  of  continuity  since  eliminating  a  control  point  can  be 
the  result  of  imposing  algebraic  equations  that  relates  the  derivatives  at  a  certain  breakpoint. 
From  the  bidirectional  structure  used  in  Eq.  1,  a  surface  segment  which  has  cubic  interpolation 
along  u  (p  =  3,  n  =  3,  r  + 1  =  8 )  and  a  linear  interpolation  along  v  ( q  =  1,  m  =  1,  5  + 1  =  4 ),  should 

have  (n  +  l)x(m  +  l)  =  8  control  points;  this  is  regardless  of  whether  the  surface  is  two-  or  three- 
dimensional.  Manipulation  of  the  B-spline  surface  of  Eq.  1  shows  that  these  eight  control  points 
are  the  result  of  using  the  alternate  basis  set  l,u,v,uv,u2,u2v,u3,u3v .  That  is,  B-spline 
representation  and  the  formulas  of  Eq.  3  do  not  allow  for  the  use  of  the  basis  set  1  ,u,v,uv,u2 ,u3 
which  can  be  effectively  used  to  develop  a  shear  deformable  beam  model.  If  a  cubic  interpolation 
is  used  for  both  u  and  v  (thin  plate),  the  B-spline  representation  will  require  16  control  points 
because  the  expansion  must  include  all  terms  mV;  k,l  =0,1, 2,3  regardless  of  whether  the  shape 
of  deformation  of  the  plate  is  simple  or  complex;  one  must  strictly  follow  the  B-spline  rigid 
structure.  This  can  be  of  disadvantage  in  the  analysis  since  such  a  geometric  representation  can 
unnecessarily  increase  the  dimension  of  the  analysis  model  and  leads  to  the  loss  of  the  flexibility 
offered  by  the  FE  method  or  modal  analysis  techniques.  As  the  degree  of  the  polynomial 
interpolation  increases,  the  problem  gets  even  worse.  Another  important  and  interesting  issue 
with  regard  to  the  use  of  B-spline  as  an  analysis  tool  is  capturing  discontinuities;  this  is  discussed 
in  the  following  section. 
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3.  STRUCTURAL  AND  NON-STRUCTURAL  DISCONTINUITIES 

As  previously  mentioned,  in  the  case  of  structural  discontinuity,  there  is  no  relative  rigid  body 
motion  at  the  discontinuity  node;  all  the  relative  displacements  are  due  to  deformations.  Example 
of  structural  discontinuity  is  at  point  C  in  Fig.l.  At  the  structural  discontinuity  node,  in  the  case 
of  the  planar  system  shown  in  Fig.  1,  there  is  only  one  state  of  strains.  Non-structural 
discontinuity,  on  the  other  hand,  allows  for  relative  rigid  body  rotation.  At  point  O  of  Fig.  1, 
there  can  be  two  different  states  of  strains  because  of  the  rigid  body  mode.  Nonetheless, 
structural  and  non-structural  discontinuities  can  be  classified  as  C° .  They  are,  however, 
fundamentally  different  from  the  analysis  point  of  view  since  they  lead  to  completely  two 
different  joint  models  that  have  different  numbers  of  degrees  of  freedom.  Interestingly,  one  of 
these  types  of  discontinuities  can  be  captured  by  B-spline  recurrence  formula,  while  the  other 
cannot  be  captured.  More  interestingly,  the  discontinuity  type  captured  by  B-spline  is  the  non- 
structural  discontinuity  that  characterizes  many  mechanical  system  applications.  Modeling 
structural  discontinuity  that  characterizes  both  mechanical  and  structural  systems  requires  the  use 
gradient  transformation  which  is  crucial  in  the  ANCF  geometric  representation. 

As  previously  mentioned,  reducing  the  knot  multiplicity  by  one  at  a  breakpoint  leads  to 
C°  continuity  and  to  the  elimination  of  one  control  point.  This  equivalent  to  the  MBS  pin  joint 
constraint  definition  in  planar  analysis  and  to  the  MBS  spherical  joint  definition  in  the  spatial 
analysis.  This  type  of  C°  continuity  that  is  captured  by  B-spline  is  of  the  non-structural 
discontinuity  type  which  leads  to  a  rigid  body  mode  and  to  a  non-unique  state  of  the  strain  at  the 
discontinuity  node.  The  B-spline  recurrence  formula  structure  leads  automatically  to  this  type  of 
discontinuity.  The  structural  discontinuity,  while  it  is  also  of  the  C°  type,  requires  additional 
algebraic  equations  in  order  to  define  a  unique  strain  state  by  eliminating  the  relative  rotation  at 
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the  joint  definition  point.  These  algebraic  equations  can  be  used  to  eliminate  other  control  points, 
and  such  elimination  is  not  embedded  in  the  rigid  B -spline  geometry  representation.  Only  one 
type  of  C°  continuity  that  can  be  achieved  by  using  B -spline  formula;  reducing  the  knot 
multiplicity  by  one  in  B-spline  representation  does  not  capture  structural  discontinuity. 

4.  GENERALITY  OF  ANCF  GEOMETRY 

While  B-spline  geometry  can  always  be  converted  to  ANCF  geometry,  the  converse  is  not  true. 
ANCF  geometry  does  not  impose  restriction  on  the  basis  functions  that  must  be  included  in  the 
interpolating  polynomials.  This  allows  for  developing  finite  elements  that  have  less  coordinates 
as  compared  to  those  developed  using  the  B-spline  representation.  Furthermore,  ANCF  geometry 
can  be  used  to  model  both  structural  and  non- structural  discontinuities  (Shabana  and  Mikkola, 
2003;  Shabana,  2010;  Hamed  et  al.,  2011),  while  the  rigid  recurrence  B-spline  representation 
cannot  be  used  to  model  structural  discontinuities  in  a  straight  forward  manner.  The  basic 
differences  between  ANCF  and  B-spline  geometries  are  demonstrated  in  this  section  using  a 
planar  beam  example.  The  displacement  field  of  the  shear  deformable  beam  used  in  this  section 
can  be  written  as  r(x,  y)  =  S(x,  y)e(t) ,  where  x  and  y  are  the  element  spatial  coordinates,  t  is 

time,  S  is  the  element  shape  function  matrix,  and  e  is  the  vector  of  the  element  nodal 
coordinates.  The  shape  function  matrix  for  the  element  considered  in  this  section  is  defined  as 

SJ  =  [,',I  s2I  s3I  ,s-4I  s5I  ,s6l]  (4) 

where  the  shape  functions  sj,i  =  1, 2, . . ., 6  are  defined  as  (Omar  and  Shabana,  2001) 
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s,=l-3<f +2<f,  !,=;(£- 2?+?),  s,=ln(  1-f),' 
s4  =  3f--2£\  ss  =  l(-f-  +^),  ss=l$ri 


In  this  equation,  £  =  x//,  7]  =  y/l.  ANCF  finite  elements  employ  gradient  vectors  as  nodal 

coordinates.  For  the  element  used  in  this  section,  the  vector  of  nodal  coordinates  is  defined  as 


e  = 


(r)i  (5r/&)'  {fcldy)\  (rf2  /dxf2  (dr/dyf2 


(6) 


where  (  )  ,  k  =  1, 2  indicates  variables  evaluated  at  node  k  of  the  element.  Note  that  the  element 

defined  by  the  preceding  equations  is  based  on  a  cubic  interpolation  for  x  and  a  linear 
interpolation  for  y  .  This  element  has  been  widely  used  in  the  analysis  of  large  deformation 
problems. 

The  finite  element  described  in  this  section  is  an  example  of  ANCF  elements  that  cannot 
be  converted  to  B -spline  representation.  This  element  is  based  on  a  polynomial  expansion  that 
does  not  have  the  two  basis  functions  x2y  and  x3y  .  These  terms  can  be  systematically  included 
in  ANCF  geometry  by  adding  nodal  coordinates  allowing  for  converting  B-spline  representation 
to  ANCF  representation.  Similar  comments  apply  to  ANCF  thin  plate  elements  that  do  not  have 
to  include  all  the  basis  functions  xkyl ;  k,l  =  0,1, 2,3  .  This  flexibility  offered  by  ANCF  geometry 
allows  for  developing  finite  elements  that  have  smaller  number  of  coordinates  compared  to  those 
elements  developed  by  B-spline  geometry. 

One  can  also  show  that  ANCF  finite  elements  can  describe  structural  and  non- structural 
discontinuities.  Non- structural  discontinuities  that  allow  for  large  rigid  body  rotations  can  be 
described  using  a  C°  model  obtained  by  imposing  constraints  on  the  position  coordinates  only. 
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For  example  if  two  elements  i  and  j  are  connected  by  pin  joint  at  a  node,  one  can  apply  the 

algebraic  equations  r'  =  rJ  at  this  node.  These  algebraic  equations  can  be  imposed  at  a 
preprocessing  stage  to  eliminate  the  dependent  variables  and  define  FE  mesh  that  has  a  constant 
mass  matrix  and  zero  Coriolis  and  centrifugal  forces  despite  the  finite  rotations  allowed  between 
the  finite  elements  of  the  mesh.  As  previously  mentioned,  non-structural  discontinuities  can  also 
be  described  using  B-spline  geometry  by  reducing  the  knot  multiplicity  at  the  joint  node  by  one. 
Note  that  in  the  case  of  non-structural  discontinuities  no  constraints  are  imposed  on  the  gradient 
vectors,  and  therefore,  the  state  of  strain  is  not  unique  at  the  joint  node.  Each  of  the  Lagrangian 

strains  =  (rjr.  -l)/2,  syy  =  (rjry  -l)  jl ,  and  sxx  =  r',ry/2  have  two  values  at  the  joint  node; 
one  defined  on  element  i  and  the  other  is  defined  on  element  j .  Here 

rv  =  dr/dx,  rv  =  dr/dy,  r,  =  dr/dz . 


The  concept  of  degrees  of  freedom  widely  used  in  mechanics  is  not  considered  in 
developing  the  recurrence  relationships  on  which  B-spline  and  NURBS  geometry  are  based.  This 
represents  another  serious  limitation  when  these  computational  geometry  methods  are  used  as 
analysis  tools;  as  evident  by  the  fact  that  B-spline  geometry  cannot  describe  structural 
discontinuities.  This  type  of  discontinuities,  while  it  remains  of  the  C°  continuity  type,  requires 
imposing  additional  constraints  on  the  gradients;  these  constraints  cannot  be  captured  by  the  B- 
spline  recurrence  formula  since  they  require  the  elimination  of  additional  vectors.  In  the  case  of 
B-spline,  C°  continuity  is  achieved  by  reducing  the  knot  multiplicity  by  one,  and  this  eliminates 
one  control  point  leading  to  the  definition  of  a  pin  joint  (non-structural  discontinuity).  ANCF 
geometry,  on  the  other  hand,  allows  for  imposing  constraints  on  the  gradients  using  the  tensor 

transformation  (<9r/<9xj)  =  (<9r/<9x2)  A  ,  where  x2  =  [x{  _yjr  and  x2  =  [x2  y2f  are  two  sets  of 
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coordinate  lines,  and  A  is  the  matrix  of  coordinate  line  transformation.  Using  this  tensor 
gradient  transformation,  the  structural  discontinuities  can  be  systematically  modeled  using 
ANCF  finite  elements  (Shabana  and  Mikkola,  2003;  Shabana,  2010).  For  example,  if  the  axis  of 
a  beam  element  j  makes  an  angle  a  with  the  axis  of  another  element  i  and  the  two  elements 
are  rigidly  connected  at  a  node,  the  structural  discontinuity  conditions  that  eliminate  all  the 
relative  rigid  body  displacements  can  be  written  at  the  joint  node  as 

dr'  dr'  drj  drj  cos  a  sin  a 

dx  dy'  dxj  dyj  _  -sin  a  cos  a 

These  six  scalar  algebraic  equations  can  be  used  to  eliminate  a  position  coordinate  vector  and 
two  gradient  vectors,  defining  a  unique  strain  state  at  the  node  of  connectivity  between  the 
ANCF  finite  elements.  The  algebraic  conditions  of  Eq.  7  that  allow  ANCF  finite  elements  to 
describe  two  types  of  C°  discontinuity  cannot  be  automatically  captured  by  the  B -spline 
recurrence  formula.  Note  that  these  algebraic  conditions  are  linear  in  the  ANCF  finite  element 
coordinates,  and  therefore,  they  can  be  applied  using  a  standard  FE  assembly  procedure  at  a 
preprocessing  stage  of  the  analysis.  These  conditions  lead  to  a  relative  motion,  between  the  finite 
elements,  that  is  pure  deformation  displacement. 

It  is  important  to  point  out  that  the  coordinate  line  transformation  of  Eq.  7  need  to  be 
applied  only  in  the  case  of  structural  discontinuities.  Such  a  transformation  is  not  required  in  the 
case  of  non-structural  discontinuities  despite  the  fact  that  the  elements  can  have  arbitrary 
orientation  relative  to  each  others.  This  is  due  to  facts  that  non-structural  discontinuity  does  not 
impose  constraints  on  the  gradient  vectors,  ANCF  geometry  is  invariant  under  an  arbitrary  rigid 
body  displacement,  and  each  ANCF  finite  element  has  its  own  independent  parameters. 
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Similarly,  in  the  case  of  curved  shapes  as  in  belt  applications,  the  coordinate  transformation  used 
in  Eq.7  is  not  required  since  there  are  no  structural  discontinuities  and  the  gradients  at  all  points 
on  the  belt  can  be  defined  with  respect  to  the  same  coordinate  lines.  In  the  curved  sections  of  the 
belt,  one  must  provide  the  appropriate  values  of  the  gradient  vectors  that  define  the  correct 
shape.  Recall  that  a  curve,  regardless  of  its  shape,  requires  only  one  parameter;  while  a  surface, 
regardless  of  its  shape,  requires  only  two  parameters. 


5.  THREE-DIMENSIONAL  NON-STRUCTURAL  DISCONTINUITIES 

In  this  section,  it  is  shown  how  fully  parameterized  ANCF  three-dimensional  finite  elements  can 
be  used  to  develop  spatial  joint  models  that  allow  large  relative  rigid  body  rotation  between  the 
finite  elements.  ANCF  finite  elements  connected  by  this  joint  can  be  assembled  using  linear 
connectivity  conditions  leading  to  FE  mesh  that  has  a  constant  mass  matrix  and  zero  Coriolis  and 
centrifugal  forces.  The  fully  parameterized  three-dimensional  ANCF  beam  element  is  used  in 
this  investigation  to  demonstrate  the  development  of  such  joint  models.  The  displacement  field 
of  the  element  can  be  written  as  r(x,  y,z)  =  S(jc,  y,z)e(t)  where  x,y,  and  z  are  the  element 

spatial  coordinates.  The  shape  function  matrix  of  this  element  is  defined  as  (Yakoub  and 
Shabana,  2001;  Shabana,  2008) 

S  =  [s,I  s2l  ,s'3I  ,s4I  ,v5I  ,s'fiI  ,s-7I  ,s8I  ]  (8) 


where  the  shape  functions  si , ,  i  =  1, 2, . . . ,  8  are  defined  as 


Sl=l-3^  +  2f,  s2=l(Z- 2?+?), 

S3  *4  =l(s-&),  s5=3%2-2%\> 

s6=l(-^2+^3),  s7=lZ;ri,  st=l$g 


(9) 
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In  this  equation,  £  =  x/l,  ij  =  y/l,  g  =  z/l .  The  element  has  two  nodes;  each  node  has  12  nodal 

coordinates  defined  by  the  vector  ek  =  |^rA/  r';'  rn  r  r/J  ,  where  k  is  the  node  number. 

The  ANCF  finite  element  defined  by  Eqs.  8  and  9  captures  the  cross  section  deformation  and  its 
coupling  with  extension  and  bending.  Therefore,  this  element  can  be  used  to  develop  general 
models  for  belt  drives  and  rubber  tracked  vehicles. 

The  three-dimensional  beam  element  presented  in  this  section  is  another  example  that  can 
be  used  to  demonstrate  the  generality  of  the  ANCF  geometry.  This  element  is  based  on  cubic 
interpolation  in  x  and  linear  interpolation  in  y  and  z  ■  Nonetheless,  one  can  show  that  the  four 

basis  functions  x2y,  x3y,  x2z,  x3z  are  not  used  in  developing  the  displacement  field  of  this 
widely  used  ANCF  beam  element.  Therefore,  the  geometry  of  this  element  cannot  be  converted 
to  B-spline  volume  geometry.  These  missing  basis  functions  can  be  systematically  included  in 
the  development  of  another  ANCF  finite  element  that  can  be  converted  to  B-spline  volume 
geometry.  However,  such  an  element  will  lead  to  50%  increase  in  the  number  of  the  element 
nodal  coordinates. 

A  planar  pin  joint  between  rigid  or  flexible  bodies  is  an  example  of  C°  continuity,  as 
previously  discussed.  A  pin  joint  between  two  rigid  bodies  in  the  spatial  analysis  also  allows  for 
only  one  degree  of  freedom,  which  is  a  relative  rotation  about  the  joint  axis.  Since  the  pin  joint 
eliminates  five  degrees  of  freedom  in  the  rigid  body  analysis,  its  formulation  requires  five 
algebraic  constraint  equations  that  eliminate  three  relative  translation  displacements  and  two 
relative  rotations  between  the  two  bodies.  In  the  case  of  flexible  bodies,  an  infinitesimal  volume 
can  have  12  modes  of  displacements;  three  rigid  body  translations,  three  rotations,  and  six 
deformation  modes.  In  this  section,  a  new  model  of  pin  joint  between  ANCF  finite  elements  is 
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introduced.  The  formulation  of  this  pin  joint  between  elements  i  and  j  employs  the  following 
six  scalar  equations  defined  at  the  joint  node: 


=  rj 


(10) 


Here  a  is  the  coordinate  line  that  defines  the  joint  axis;  a  can  be  x,  y ,  or  z  or  any  other 
coordinate  line  as  discussed  later  in  this  section.  The  six  scalar  equations  of  Eq.  10  eliminate  six 
degrees  of  freedom;  three  translations,  two  rotations,  and  one  deformation  mode.  This  joint 
model  ensures  C1  continuity  with  respect  to  the  coordinate  line  a  and  C°  continuity  with 
respect  to  the  other  two  parameters.  It  follows  that  the  Lagrangian  strain  component 

8 aa  =(rJrff  — 1)/2  is  continuous  at  the  joint  definition  point,  while  the  other  five  strain 
components  can  be  discontinuous. 

While  the  algebraic  constraint  equations  of  a  pin  joint  between  two  rigid  bodies  are 
highly  nonlinear.  The  algebraic  constraint  equations  of  Eq.  10  are  linear.  Therefore,  these 
equations  can  be  applied  at  a  preprocessing  stage  to  systematically  eliminate  the  dependent 
variables.  Using  these  equations,  one  can  develop  a  new  kinematically  linear  FE  mesh  for 
flexible-link  chains  in  which  the  links  can  have  arbitrarily  large  relative  rotations.  The  use  of  this 
pin  joint  model  with  ANCF  finite  elements  leads  to  a  constant  mass  matrix  and  zero  Coriolis  and 
centrifugal  forces. 

As  previously  mentioned  in  this  paper,  B-spline  geometry  can  describe  the  type  of  non- 
structural  discontinuity  discussed  in  this  section.  Nonetheless,  if  an  arbitrary  axis  of  a  pin  joint  is 
to  be  used  in  the  analysis,  the  use  of  B-spline  geometry  can  be  difficult  given  the  rigid  structure 
of  the  B-spline  recurrence  formula.  In  order  to  be  able  to  choose  an  arbitrary  axis  of  rotation  for 
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the  pin  joint,  one  must  be  able  to  define  the  gradient  vector  in  the  direction  of  a  coordinate  line 
along  this  axis  of  rotation.  Such  a  definition  can  be  easily  made  using  ANCF  geometry  using  the 
gradient  tensor  transformation.  Let  u,v  ,  and  w  be  another  set  of  parameters;  one  of  them  can  be 

used  to  define  the  joint  axis.  It  follows  that  [rf  r  rw]  =  [rv  rv  r  ]  A ,  where  A  is  the 

constant  matrix  of  coordinate  transformation  defined  as 


dx 

dx 

dx 

du 

dv 

dw 

dy 

dy 

dy 

du 

dv 

dw 

dz 

dz 

dz 

du 

dv 

dw 

(11) 


The  fact  that  this  matrix  is  constant  allows  having  linear  pin  joint  connectivity  conditions  when 
ANCF  finite  elements  are  used  (Shabana  and  Mikkola,  2003;  Shabana,  2010). 


6.  NUMERICAL  EXAMPLES 

In  this  section,  three-dimensional  belt  drive  and  chain  examples  are  used  to  demonstrate  the 
implementation  of  the  concepts  discussed  in  this  paper.  The  degree  of  continuity  at  the  element 
interfaces  can  be  applied  at  a  preprocessing  stage  in  order  to  eliminate  the  dependent  variables 
leading  to  a  finite  element  mesh  that  has  a  constant  mass  matrix  and  zero  Coriolis  and  centrifugal 
forces.  The  fully  parameterized  three-dimensional  ANCF  beam  element  discussed  in  this  paper  is 
used  in  modeling  the  belt  and  the  chain  used  in  this  section.  The  following  three  different  models 
are  considered  in  this  section: 

1.  A  finite  element  belt  model  with  C1  continuity.  The  geometry  of  this  model  can  be 
defined  in  the  initial  configuration  using  a  smooth  curve  that  defines  the  centerline  of  the 
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belt.  This  model,  referred  to  in  this  section  as  the  C1  belt  model,  leads  to  continuous 
gradients  and  strains  at  the  element  interfaces. 

2.  A  finite  element  belt  model  that  ensures  continuity  of  r  (C1  continuity  along  y  ),  but  it 

has  C°  continuity  along  x  and  z  .  In  the  initial  configuration,  the  centerline  of  the  belt  is 
continuous.  This  model  is  referred  to  as  the  C°/Cl  belt  model. 

3.  A  finite  element  flexible-link  chain  model  that  ensures  continuity  of  r  (C1  continuity 


along  y  ),  but  it  has  C°  continuity  along  x  and  z  .  This  model  differs  from  the  previous 
two  models  due  to  the  fact  that  the  chain  centerline  is  not  continuous  at  the  initial 
configuration  because  rx  and  r  are  not  continuous.  This  model  will  be  referred  to  in  this 
section  as  the  chain  model. 

The  belt  and  chain  drive  mechanisms  considered  in  this  section  are  assumed  to  consist  of  two 
pulleys  and  a  flexible  belt  or  a  chain  as  shown  in  Figs.  2  and  3.  Figure  3  shows  the  discontinuity 
of  the  gradients  at  the  joints  in  the  initial  configuration.  The  pulleys  are  connected  to  the  ground 
using  revolute  joints.  The  span  length  in  the  reference  configuration  is  assumed  to  be  0.12  m  in 
the  case  of  the  belt  model  and  0.128  m  in  the  case  of  the  chain  model.  In  the  examples 
considered  in  this  section,  the  angular  velocity  of  the  driving  pulley  is  specified  by  the  following 
equation: 


0 

t  <  0.3 

CO,  =< 

100«-03) 

0.3<t  <1.0 

1 

0.7 

100 

f  >1.0 

(12) 


where  a\  is  expressed  in  rad/s,  and  t  is  the  time  expressed  in  seconds.  The  mass  moment  of 

2 

inertia  of  the  driven  pulley  about  its  axis  of  rotation  is  0.0075  kg/m".  In  order  to  introduce 
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tension  in  the  system,  the  driving  pulley  is  displaced  in  the  X  direction  until  the  span  length 
reaches  0.18  m  for  the  three  models.  In  order  to  avoid  oscillations  at  the  beginning  of  the 
simulation,  this  pulley  displacement  is  achieved  over  a  period  of  0.2  s.  A  resistance  moment 
defined  by  the  following  equation  is  also  applied  to  the  driven  pulley: 

0  t  <  0.3 

M  =  j-l.l(r~Q'3)  0.3  <  t  <1.0  (13) 

0.7 

-1.1  r  >1.0 

The  two  pulleys  are  assumed  to  have  the  same  radius  and  width  of  0.058  and  0.01  m, 
respectively.  A  compliant  force  model  is  used  to  describe  the  belt/pulley  interaction.  The 
stiffness  and  damping  coefficients  used  in  the  belt/pulley  contact  force  model  are  given 

respectively  by  9x10  N/m  and  2  x  1 0  N.s/m  .  Tangential  friction  forces  are  also  introduced 
using  a  coefficient  of  dry  friction  of  1.2  (Leamy  and  Wasfy,  2002).  The  friction  parameter  that 
defines  the  slope  in  the  transition  region  is  assumed  to  be  106  N.s/m3  (Dufva  et  al,  2007).  The 
belt  is  modeled  using  20  ANCF  three-dimensional  beam  elements.  Incompressible  Neo- 
Hookean  constitutive  model  with  nonlinear  damping  model  is  used  to  model  the  flexible  belt. 
The  belt  is  assumed  to  have  a  rectangular  cross-section  of  dimensions  0.01  x  0.004  m  and  density 
of  3500  kg/m  .  The  incompressible  Neo-Hookean  model  constant  is  assumed  to  be 

jus  =2xl06N/m2,  the  incompressibility  constant  is  assumed  to  be  k  =  10*  N/m2,  and  the 

dilatation  and  deviatoric  dissipation  factors  used  for  the  damping  model  are  assumed  to  be  10  4 
and  5  x  10  5 ,  respectively. 

Figure  4  shows  the  angular  velocity  of  the  driving  and  driven  pulleys  for  both  the  C1 , 
C°/Cl  ,  and  chain  models.  The  results  presented  in  this  figure  show  that  the  angular  velocity  of 
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the  driven  pulley  in  the  case  of  the  C1  belt  model  is  higher  than  that  of  C°/C'  belt  and  chain 
models.  This  can  be  attributed  to  the  fact  that  in  the  case  of  the  C0/ C1  and  chain  models,  some 
of  the  loads  such  as  bending  moment  about  Y  axis  are  not  transferred  between  elements.  Figures 
5-7  show  the  configurations  of  the  belt  centerline  for  the  three  models  at  time  1  s.  These  figures 
show  that  the  gradients  rv  and  r  are  discontinuous  in  case  of  the  C0/ C1  belt  and  chain  models 

and  they  are  continuous  in  case  of  the  C1  belt  model. 

Figure  8  shows  a  measure  of  the  cross  section  deformation  along  the  belt  centerline  at 
time  1.9s.  This  measure  is  defined  by  Nanson’s  formula  that  can  be  used  to  calculate  the  ratio 
between  the  areas  in  the  current  and  reference  configurations.  Nanson’s  formula  is  defined  as 
(Ogden,  1984;  Shabana,  2008) 


where  a  and  A  are,  respectively.  The  area  in  the  current  and  the  reference  configurations,  J  is 
the  determinant  of  the  matrix  of  position  vector  gradients  J  ,  and  nn  is  the  unit  vector  normal  to 

the  area.  The  results  of  Fig.  8  show  that  the  area  ratio  is  continuous  in  the  case  of  the  C1  belt 
model  while  it  is  not  continuous  in  case  of  the  C°/C'  and  chain  models.  Figures  9-11  show 
comparison  of  the  normal  strains  at  the  element  interface  points  along  the  belt  centerline  for  the 
three  models.  The  results  obtained  in  this  investigation  show  that  all  the  strain  components, 
which  are  functions  of  the  gradient  vectors,  are  continuous  in  the  case  of  the  C1  belt  model; 
while  in  the  case  of  the  C°/Cl  and  chain  models  only  s  is  continuous.  Figure  12  shows  the 
Lagrangian  shear  strain  component  sxz  along  the  belt  centerline.  The  results  of  this  figure  show 
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sx_  to  be  discontinuous  in  the  case  of  the  C°/Cl  and  chain  models  because  it  is  function  of  both 
the  gradient  vectors  rt  and  r, ,  which  are  discontinuous  when  these  models  are  used. 

7.  SUMMARY  AND  CONCLUSIONS 

This  paper  addresses  the  important  issue  of  using  computational  geometry  methods  such  as  B- 
spline  and  NURBS  as  analysis  tools.  B-spline  and  NURBS  employ  recurrence  formulas  that 
allow  changing  the  degree  of  continuity  at  a  breakpoint  by  adjusting  the  knot  multiplicity  at  this 
point.  As  demonstrated  in  this  paper,  the  recurrence  formula  has  several  drawbacks  when  B- 
spline  representation  is  used  as  an  analysis  tool.  Because  the  recurrence  formula  does  not  provide 
flexibility  for  choosing  the  basis  functions,  B-spline  representation  can  lead  to  significantly 
larger  number  of  coordinates  and  a  higher  dimensional  model.  This  fact  was  used  to  demonstrate 
the  generality  of  the  ANCF  geometry.  While  B-spline  geometry  can  always  be  converted  to 
ANCF  geometry,  the  converse  is  not  true.  It  is  also  shown  that  the  B-spline  recurrence  formula 
cannot  be  used  to  model  structural  discontinuity  in  a  straight  forward  manner.  While  structural 
discontinuities  are  of  the  C°  type,  they  cannot  be  captured  in  the  B-spline  representation  by 
reducing  the  knot  multiplicity  by  one.  This  reduction  of  the  knot  multiplicity  is  equivalent  to 
elimination  of  the  relative  translation  only;  and  such  a  reduction  in  the  knot  multiplicity  leads  to 
a  rigid  body  mode  that  defines  the  conditions  of  a  pin  joint.  In  the  case  of  structural 
discontinuities,  on  the  other  hand,  the  C°  B-spline  representation  does  not  eliminate  the  rigid 
body  mode  since  additional  algebraic  constraint  equations  are  required  in  order  to  eliminate  the 
relative  rotations  between  two  segments.  The  paper  also  presents  a  new  three-dimensional  pin 
joint  model  that  leads  to  linear  connectivity  conditions  and  constant  mass  matrix  when  used  with 
ANCF  finite  elements.  The  implementation  of  this  new  model  is  demonstrated  using  a  flexible  - 
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link  chain  example.  The  limitations  identified  in  this  paper  when  B-spline  geometry  is  used  as 
analysis  tool  suggest  the  use  of  the  I-CAD-A  approach  in  which  a  constant  transformation  can  be 
developed  to  convert  CAD  geometry  to  FE  mesh.  It  should  be  also  clear  that  NURBS  geometry 
has  the  same  limitations  as  B-spline  representation  when  used  as  an  analysis  tool. 
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Figure  1  Structural  and  non-structural  discontinuities 
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Figure  2  Initial  configuration  of  the  belt  drive  mechanism  for  both  CH/ C1  and  C1  models 
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Figure  4  Angular  velocity  of  the  driving  and  driven  pulleys 
(— & — Driving  pulley,— •—  C1  belt  (driven  pulley),  C°/Cl  belt  (driven  pulley), 

-  *  -  chain  (driven  pulley)) 
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Figure  5  Centerline  of  the  C1  belt  model  at  time  Is 


Figure  6  Centerline  of  the  C°/C'  belt  model  at  time  Is 
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Figure  9  Axial  strain  sxx  along  the  belt  centerline  at  time  1.9s 
(  — •—  C1  belt,  A  CVC1  belt,  chain) 


Figure  10  Normal  strain  e  along  the  belt  centerline  at  time  1.9s 
(  — •—  C1  belt,  A  C0/ C1  belt, chain) 


